Effect of Particle Shape and Charge on Bulk Rheology of 

Nanoparticle Suspensions 



David R. Heine 

Science and Technology Division, Corning Incorporated, 
SP-TD-01-1, Corning, New York 14831, USA 



Matt K. Peterseqj and Gary S. Grest 
Sandia National Laboratories, Albuquerque, New Mexico 87185, USA 

(Dated: April 15, 2010) 

Abstract 

The rheology of nanoparticle suspensions for nanoparticles of various shapes with equal mass 
is studied using molecular dynamics simulations. The equilibrium structure and the response to 
imposed shear are analyzed for suspensions of spheres, rods, plates, and jacks in an explicit solvent 
for both charged and uncharged nanoparticles. For the volume fraction studied, v f = 0.075, the 
uncharged systems are all in their isotropic phase and the viscosity is only weakly dependent on 
shape for spheres, rods, and plate whereas for the jacks the viscosity is an order of magnitude 
larger than for the other three shapes. The introduction of charge increases the viscosity for all 
four nanoparticle shapes with the increase being the largest for rods and plates. The presence 
of a repulsive charge between the particles decreases the amount of stress reduction that can be 
achieved by particle reorientation. 
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I. INTRODUCTION 

With increasing capabilities to synthesize nanoparticles with a wide variety in composi- 
tions, shapes and sizes there has been an increasing interest in how these properties influence 
the behavior of the bulk material. Understanding what nanoparticle properties most strongly 
influence the bulk behavior and to what extent the bulk properties vary can provide a signif- 
icant advantage in a number of industries that process powders, suspensions, and solutions 
containing nanoscale particles. One such property is the shape of the nanoparticle which 
even in naturally occurring materials can vary drastically from elongated rods to near perfect 
spheres or amorphous blobs. To what extent do these various shapes impact the rheology 
of nanoparticle suspensions? How much of a difference in nanoparticle shape is required 
to use tailored nanoparticle shapes as viscosity modifiers? We address these questions by 
simulating the rheological behavior of suspensions of nanoparticles with different shapes. 

Multiple approaches have been used to study the rheology of suspensions including Brown- 
ian dynamics,- 1 ^ dissipative particle dynamics,- - - and stochastic rotation dynamics (SRD).-^- 
Almost all of these approaches treat the nanoparticles as perfectly spherical particles in ei- 
ther an implicit solvent or an explicit solvent consisting of spherical particles. For SRD, 
the implicit solvent particles are treated as point masses. Some work has been done us- 
ing quadrupoles to describe suspensions of platelets representing clay particles instead of 
spheres^ - — but this approach does not allow for a wide range of particle shapes. Others 
have described particles of arbitrary shapes using composites of sub-particles held together 
as rigidii^ or flexible bodies^. This is the approach we use here to study the effect of 
nanoparticle shape on suspension rheology for the four nanoparticle shapes shown in Fig. 
[TJ In addition, we also consider how the introduction of charge on the nanoparticle affects 
their suspension rheology in the presence of both low and high salt solutions. 

In Section [Til we define the systems studied including the interaction potentials between 
nanoparticles and the system properties. We also describe the techniques used to simulate 
the equilibrium and sheared systems. In section [TTT1 we present our results for the equilibrium 
structure and diffusion coefficient for nanoparticle suspensions for the four nanoparticle 
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FIG. 1: Snapshots of individual particles for (a) jacks consisting of three 25a long and 1.2a wide 
arms, (b) rods that are 22.5a long and 2a wide, (c) plates with each side 12a long, and (d) spheres 
with diameters of 5.6a. 

shapes. We determine the relationship between the diffusion coefficient and ionic strength 
for charged spheres and jacks as we approach the no salt limit. Results for the orientational 
order parameter for rods and plates for different ionic strength are also presented. In section 
IIVI we present results for the viscosity of the sheared suspensions. A brief summary of our 
results and conclusions are presented in Section [V] 

II. MODEL DETAILS 

Nanoparticles in the shape of rods, plates, jacks and spheres were each modeled as rigid 
composites of 135 interaction sites each of mass m as shown in Fig. [TJ The uncharged double 
wide jacks discussed here as well as extended jacks consisting of single row of interaction 
sites have previously been studied by Petersen et al.— The background solvent consists 
of Lennard- Jones particles of mass m. The interaction potential for both the solvent and 
nanoparticle sites is the Lennard- Jones (LJ) 6-12 potential, 
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where e and a are the LJ units of energy and length and are the same for both the nanopar- 
ticle interaction sites and the solvent. The solvent-solvent and solvent-nanoparticle site 
interactions were cutoff at 3a while the nanoparticle-nanoparticle site interactions were cut- 
off at 2 1 / 6 cr. The shorter cutoff for the nanoparticle-nanoparticle site interaction results in a 
purely repulsive, short range interaction between nanoparticles, thereby promoting nanopar- 
ticle dispersion. The equations of motion were integrated with a velocity- Verlet algorithm 
with a time step At = 0.005r, where r = a(m/e) 1 ^ 2 . The density of the nanoparticles 
was nominally the same as that of the background solvent to avoid the consideration of 
sedimentation. The rigid-body motion was determined by the forces on the nanoparticle's 
composite sites. These forces produced center-of-mass translation and rotation about the 
nanoparticle center-of-mass. All of the simulations for the uncharged systems were carried 
out at a nanoparticle volume fraction of v f = 0.075. This volume fraction was chosen since 
it was large enough that the suspension viscosity r\ was measurably larger than that of the 
solvent viscosity, r] s = 1.01 ± 0.03m/r<7,— and low enough that at least for the uncharged 
systems the nanoparticles were in their isotropic phase. Although each shape is composed 
of the same number of sub-particles, the surface area and hydrodynamic radius varies for 
each type, thus giving rise to shape-dependent rheological behavior. 

The effect of introducing charge on the nanoparticles was also studied by adding a DLVO 
interaction between all interactions sites on the nanoparticles, 



where the prefactor A is set to l.Oe and the cutoff is at r c = 10a. We compare high screening 
and low screening to the no charge case for all four systems by setting the inverse Debye 
length to k — 0.34a -1 for the high screening case and k = 0.086a -1 for the low screening 
case. These values correspond to Debye lengths of d/2 and 2d, respectively, where d = 5.8a 
is the diameter of the spherical composite particle. Upon adding the charges to each system, 
we re-equilibrate the system volume using a Nose/Hoover barostat which leads to a slight 
increase in the system volume. On average, the volume increases by 1% for k = 0.34a _1 and 
by 3% for n = 0.086a -1 . For the spheres and jacks, we also simulated larger values of k to 
explore the crossover to the limit of no salt (k, — > oo). It should be noted that although the 
DLVO potential provides an accurate representation of the forces between charged particles 
at large separations, the theory does not take into account the distortion of the electrical 
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System L/a N solv Nj 



part 



rods 



71.96 236367 143 



jacks 72.34 236367 143 



plates 41.34 44546 27 



spheres 41.38 44550 27 



TABLE I: Simulation box length L and number of solvent particles (N so i v ) and number of nanopar- 
ticles (N part ) for the four nanoparticle shapes 

double-layer that becomes significant when a pair of particles come within two Debye lengths 
of each other. We don't expect this to be a severe limitation for our dilute solutions, but it 
will contribute some quantitative error when comparing these results to that for real particle 
suspensions. 

Given the elongation of the jacks and rods, these two systems require a larger simulation 
cell L than for the spheres and plates to avoid any issues related to the finite simulation 
size. The system sizes and number of particles in each system are summarized in Table 
[H Each simulation uses periodic boundary conditions in all directions. The temperature is 
maintained at T = e/fce using a Langevin thermostat with damping constant T = O.Olr -1 . 
For shear simulations, the thermostat is applied only in the direction perpendicular to the 
bulk flow and shear directions. 

Nonequilibrium molecular dynamics (NEMD) simulations are performed using the reverse 
perturbation method of Mulier-Plathe.— In this approach, particle momentum is exchanged 
between particles at the boundary and in the middle layer of the simulation domain to induce 
a shear velocity profile in the system. The imposed momentum flux, j z {p x ), and the velocity 
gradient, dv x /dz, are recorded every 1000 simulation steps. Using these, the viscosity for a 
given shear rate is calculated via 



The simulations are performed using the LAMMPS 16 simulation package on quad-core 2.2 
GHz AMD processors. A total of 512 cores were used for the larger jack and rod simulations 
requiring about 3 hours per million simulation steps. The sphere and plate simulations 
used primarily 64 cores which required about 4.5 hours per million simulation steps. The 
simulations were run for up to 30 million steps or 150,000 r. 



JziPx) 



dv x 
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dz 
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FIG. 2: Snapshots of the uncharged particle suspensions for (a) jacks, (b) rods, (c) plates, and (d) 
spheres at <p v f = 0.075. Solvent particles are not shown and solid particles are colored individually 
for visual clarity. 

III. STRUCTURE AND DIFFUSION RESULTS 

Snapshots of equilibrated configurations of the four uncharged suspensions are shown in 
Fig. [2j From these, we see that the jacks have a large degree of particle interactions. The 
spheres have the least amount of particle interaction due to their compact geometry. The rod 
and plates show a degree of particle alignment even in the absence of an external shearing 
force which helps to minimize the frequency of particle collisions. Snapshots of equilibrated 
configurations for the charged particles in a low salt environment are shown in Fig. |3j 
In this case the configurations of the jacks and spheres are similar to the corresponding 
uncharged systems whereas the rods tend to separate to maximize the distance between 
individual charge sites. The plates show the greatest difference in structure. The repulsive 
charge between the plates pushes apart the layered structure found in the uncharged system 
causing the plates to adopt a more randomly oriented structure. 

The change in suspension structure as a function of the shape and nanoparticle charge 
can be seen in the radial distribution function g(r). Results for g(r) measured from the 
center of each nanoparticle is shown in Fig. H] for the four nanoparticle shapes. As the salt 
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FIG. 3: Snapshots of charged particle suspensions for the low salt concentration (k = 0.086<t _1 ) 
for (a) jacks, (b) rods, (c) plates, and (d) spheres at v f = 0.073. Solvent particles are not shown 
and solid particles are colored individually for visual clarity. 

concentration decreases, the amount of local order for the jacks increases. The rods, which 
show a noticeable increase in orientation upon adding charge, show almost no dependence 
on salt concentration. The three distinct peaks in g(r) for the rods indicate a highly ordered 
structure for both low and high ionic strength. The plates show slightly higher peaks when 
reducing the salt concentration, but not as much as the spherical particle system. For the 
spherical particles, the stronger charge interaction at low ionic strength results in a larger 
primary peak indicating more local ordering than at high ionic strength. 

The change in nanoparticle shape has a strong effect on the diffusion constant D as seen 
from the results listed in Table HI1 The diffusion constant was determined from the slope 
of the mean square displacements which are shown in Fig. [5] for the two charged systems. 
Not surprisingly, the spheres show the highest diffusion constant followed by the plates, 
rods, and jacks. Increasing the salt concentration results in a 50% increase in D for the 
jacks and a twofold increase in D for the other three systems. The low values of D for the 
non-spherical nanoparticles is a consequence of their elongated structures and their inability 
to rotate freely without colliding into nearby nanoparticles. This effect is increased as the 
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FIG. 4: Radial distribution function g(r) for k = 0.34cr -1 (solid line), and k = 0.086cr _1 (dashed 
line) systems for (a) jacks, (b) rods, (c) plates and (d) spheres. 

salt concentration is reduced since the increased electrostatic repulsion leads to a larger 
effective nanoparticle size. Although D is roughly four times higher for the uncharged rod 
and plate systems compared to D for the k = 0.086a _1 systems, the values of D for the 
uncharged spheres are nearly seven times higher than that for the k = 0.086a _1 system. To 
better characterize the dependence of particle diffusivity on the solution ionic strength, we 
performed a series of simulations varying k and extracted D. The results, shown in Figure 
El indicate that D decreases fairly linearly as a function of the Debye length up to = 3cr 
suggesting that the interparticle force is sensitive to the salt concentration for this full range 
of These curves will ultimately reach a plateau at higher k' 1 as the electrostatic 

screening becomes weak enough for the particle charge to dominate the interparticle force. 
In such a case, the accuracy of the idealized screened Coulomb interaction diminishes as we 
do not take into account effects associated with the distortions of the counterion clouds. To 
further explore the role of particle shape and charge, we consider the orientational order of 
the rod and plate particles. 

The ordering in the rod and plate systems is characterized by the orientational order 
parameter, S, defined in Eq. H] where 9 is the angle between the principle axis of the particle 
and the director where the director is the preferred orientation. The brackets denote an 
average over all particles %. 

S = (3cos 2 ^ - l)/2 (4) 
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TABLE II: Diffusion constant D in units of <j 2 /t for each of the four nanoparticle shapes. 



shape k 


= 0.086O-" 1 


k = 0.34CT- 1 


uncharged 


jack 


8.9e-5 


0.00013 


0.0013 


rod 


0.00019 


0.00033 


0.00079 


plate 


0.00038 


0.00067 


0.0016 


sphere 


0.0018 


0.0041 


0.0125 




time [t] 



FIG. 5: Mean square displacement for k = 0.34<7 1 (solid line) and k = 0.086cr 1 (dashed line) 
charged systems. The colors correspond to jacks (black), rods (red), plates (green), and spheres 
(blue). 
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FIG. 6: Diffusion constant as a function of the inverse Debye screening length k for the suspension 
of jacks (black circles) and spheres (red squares). Results for the jacks are multiplied by 10 for 
clarity. 
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For each system, the orientational order parameter reached a steady state value within the 
equilibration period. The addition of charge results in an increase in S from 0.36 ± 0.10 for 
the uncharged rods to S = 0.67±0.01 for k = 0.34a- 1 and S = 0.70±0.01 for « = 0.086a" 1 . 
Whereas the addition of charge results in a decrease in S for the plates from 0.24 ± 0.12 
to S = 0.11 ± 0.05 for k = 0.34a" 1 and S = 0.16 ± 0.05 for « = 0.086a- 1 . The rods 
adopt a significantly more ordered conformation when charge is applied, but the plates do 
not. When the plates are charged, the closely stacked aggregates shown in Fig [2b are forced 
apart into a more homogeneous but disordered configuration resulting in a reduction in the 
order parameter. 

IV. RHEOLOGY RESULTS 

Each of the fully equilibrated suspensions are sheared using the Miiller-Plathe reverse 
perturbation method for shear rates ranging from 7 * 10~ 6 to 0.005 r . Snapshots of the 
various systems during shear are shown in Fig. [7] at the specified shear rates. The particle 
orientations are very comparable to the unsheared systems shown in Fig. [3] owing primarily 
to the fact that the particle charge causes the alignment of the rod particles and the breakup 
of plate aggregates even without shear. 

The momentum flux as a function of shear rate for each system simulated is shown in 
Fig. [HI In practice, we repeat the shear simulations for each system while systematically 
reducing the shear rate until the ratio of momentum flux to shear rate becomes linear. This 
indicates that we have reached the zero shear viscosity for each system. The zero shear rate 
viscosities are summarized in Table IIHI The uncharged rods, and plates have nearly equal 
viscosities that are slightly greater than that for the spheres, similar to the observations 
of Knauert et ai— from simulations at a volume fraction of <\> v f = 0.05. For each shear 
rate, the viscosity of the jacks is higher than the other shapes. Like the spheres, the jacks 
show little dependence on the ionic strength of the solution, but they do show a significantly 
higher viscosity for charged versus uncharged particles. This increase in viscosity for jacks 
compared to the other three shapes is most likely due to the very large effective volume 
occupied by the jacks, which causes them to become highly entangled as seen in Fig. 2. The 
spheres also show minimal sensitivity to the introduction of charge due to their inability to 
reduce particle interactions by orienting in the direction of flow. The small radius of gyration 
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FIG. 7: Simulation snapshots for the charged, low ionic strength (k = 0.086cr _1 ) systems at (p v f = 
0.075 under their highest simulated shear rates: (a) 7 = 0.0046r _1 for jacks, (b) 7 = 0.0020r _1 for 
rods, (c) 7 = 0.00076t -1 for plates, and (d) 7 = 0.0036r _1 for spheres. In each case, the system 
is subject to simple shear with flow in the horizontal direction and the gradient in the vertical 
direction. Solvent particles are not shown and solid particles are colored individually for visual 
clarity. 

of the spheres compared to the jacks allows them to flow past each other more readily. This 
results in a low viscosity for the suspensions of spheres relative to the other particle shapes. 
Both the rods and plates show a much higher viscosity in the low ionic strength solvent 
where the charge interaction is stronger than the high ionic strength solvent. The stronger 
charge interaction increases the effective diameter of the particles, which has a similar effect 
to increasing the volume fraction of the already constrained particles. 

We also compare the viscosities of the uncharged systems to viscosities calculated using 
a path integral technique^ in Table IIII1 The path integral approach treats the particles as 
composites of spheres as they are defined in Figure [Hand determines the intrinsic viscosity by 
first calculating the polarizability tensor of a perfect conductor of the same shape and then 
relating the intrinsic viscosity to the intrinsic conductivity for an infinitely dilute solution. 
Assuming a linear dependence of the viscosity on the particle concentration, we find that 
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FIG. 8: Log plot of momentum transferred for suspensions of various nanoparticle shapes. The 
charged suspensions have a dielectric solvent with an inverse Debye screening length of k = 0.34o" _1 
(solid lines, filled symbols) or k = 0.086a- 1 (dashed lines, open symbols) for jacks (black), rods 
(red), plates (green), and spheres (blue). 

TABLE III: Zero shear rate viscosity rj/r] s , where rj s = l.Olm/rcr is the viscosity of the neat 
solvent along with results from path integral calculations of the uncharged systems. The value in 
parenthesis gives the uncertainty in the last displayed digit. 

shape k = 0.086er -1 k = 0.34c"- 1 uncharged path integral 



jack 


31(1) 


32.8(8) 


9.8(5) 


5.1(3) 


rod 


14.9(6) 


8.2(5) 


1.7(2) 


2.5(1) 


plate 


13.6(3) 


5.7(3) 


1.9(1) 


2.2(7) 


sphere 


2.2(2) 


1.8(2) 


1.2(1) 


1.3(2) 



the agreement between the continuum approach and our simulation approach ranges from 
very good for the spherical particles to within 50% for the jacks. This deviation for the jacks 
may be due to the fact that the viscosity increased non-linearly with concentration. Since 
the effective volume of the jacks is large compared to the spheres, it is similar to being at a 
higher volume fraction such that rj > [77] c where [rf\ is the intrinsic viscosity obtained from 
the path integral approach and c is the concentration. 

The influence of shear flow on the orientation of these systems is shown in Figure |9] where 
we plot S as a function of shear rate for the charged plates and rods. In general, shear 
does not induce ordering of the particle suspensions over this range of shear rates except for 
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FIG. 9: Orientational order parameter as a function of applied shear rate for charged rods with 
k = 0.34c- 1 (black) and k = 0.086cx _1 (red) and for charged plates with k = 0.34<7 -1 (green) and 
k = 0.086a- 1 (blue). 

the highly screened plates at shear rates above 7 = O.OOlr -1 and a slight increase for the 
highly screened rods. Others^ have shown that NEMD simulations of colloidal suspensions 
can produce shear induced ordering under steady shear and that the effect is dependent on 
the simulated system size, indicating that it is an artefact of the simulation and not a real 
phenomenon. However, the shear rates shown in Figure M correspond to Deborah numbers 
ranging from De = 0.22 to De = 174. In this range, artificial shear induced ordering is found 
to occur only at reduced temperatures below about T* = 0.024 where T* = 6k B T/A. 18 The 
reduced temperature for each of our simulations is T* = 6, so we expect that the ordering 
shown in Figure [9] is not an artefact of the NEMD simulation approach. 

V. CONCLUSIONS 

Molecular dynamics simulations were carried out to study the rheology of nanoparticle 
suspensions for particle shapes of jacks, rods, plates, and spheres both charged and un- 
charged. We analyze the equilibrium structure of each system and find that the rod and 
plate systems show noticeable particle alignment which helps to minimize the frequency of 
particle collisions. We expect that even more ordering in the rod and plate systems would 
occur at higher concentrations. In addition, higher concentrations of jacks have exhibited en- 
hanced viscosity^ previous studies of spheres have found crystallization at volume fractions 
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starting around <p = O.5.— The particle diffusivity is significantly higher for the spherical 
particles than for all other particle shapes for both charged and uncharged particles. We 
also consider the response of each particle type to imposed shear. The viscosities of the 
uncharged systems are most sensitive to the ability of the particles to align under shear 
with the jack system having a significantly higher viscosity than the other particle shapes. 
Charge effects are also most prominent for shapes that can align under shear. Here, the 
presence of a repulsive charge between the particles reduces the amount of stress reduction 
that can be achieved by particle reorientation. This work characterizes the role that shape 
plays on the rheology of charged and uncharged nanoparticle suspensions. 
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